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DESCRIPTIVE  SUMMARY 


The  purpose  of  this  research  was  to  develop  a  computer  code  for  pre¬ 
dicting  the  performance  of  the  flow  in  a  dump  combustor.  It  was  felt 
that  current  elliptic  Navier-Stokes  calculations(for  example,  the 
STARPIC  program)  represent  an  inefficient  approach  to  such  problems. 
This  is  because  elliptic  calculations  are  time-consuming  and  only  a 
portion  of  the  flow  field  actually  requires  the  elliptic  approach. 
Accordingly,  a  combined  elliptic/ parabolic  scheme  was  developed  in 
which  STARPIC  was  only  used  for  the  recirculation  region.  Downstream  of 
the  recirculation  region  the  parabolized  form  of  the  Navier-Stokes 
equation  was  used.  The  code  written  for  this  region  used  an  architec¬ 
ture  similar  to  that  found  in  STARPIC.  These  programs  were  then  com¬ 
bined  producing  what  we  call  a  hybrid  code.  Calculations  were  run  in  a 
dump  combustor  having  an  area  ratio  of  approximately  four.  Comparisons 
were  made  of  the  radial  profiles  of  velocity  and  turbulence  kinetic 
energy  with  both  the  totally  elliptic  calculations  and  with  experiment. 
In  all  cases  the  hybrid  method  agreed  with  the  experiments  at  least  as 
well  as  the  elliptic  calculations,  and  in  most  cases  showed  even  better 
agreement.  Results  indicate  that  an  order  of  magnitude  savings  in  both 
computational  time  and  machine  storage  can  be  achieved  if  the  hybrid 


method  is  used. 


PUBLICATIONS,  PERSONNEL,  AND  INTERACTIONS 


At  the  present  time,  this  work  has  not  resulted  in  any  journal  publica¬ 
tions.  It  is  planned  to  publish  this  work  in  the  next  twelve  months, 
however,  most  likely  in  the  form  of  an  AIAA  conference  paper. 

Mr.  A.  A.  Hale  assisted  with  this  research,  and  with  the  financial  sup¬ 
port  provided  by  this  grant,  received  his  Master’s  Degree  in  Mechanical 
Engineering  from  the  University  in  June  of  1984.  The  title  of  his 
thesis  was  "A  New  Method  for  Calculating  the  Flow  Field  in  a  Dump  Com¬ 
bustor"  . 

On  Tuesday,  June  12,  1984  I  had  a  discussion  of  thi3  work  with  Dr.  F.D. 
Stull,  and  Dr.  R.R.  Craig  of  AFWAL/PORT  and  Dr.  P.  Vanka  of  the  Argonne 
National  Laboratory.  This  discussion  took  place  during  the 
AIAA/SAE/ASME  Joint  Propulsion  Conference  in  Cincinnati,  Ohio.  I  gave 
a  brief  review  of  the  research  and  in  the  course  of  that  discussion 
distributed  copies  of  Mr.  Hale's  thesis. 

As  a  follow-on  to  this  research  I  will  shortly  present  to  AFOSR  a  pre¬ 
proposal  for  base  flow  combustion  research  which  makes  use  of  the  com¬ 
putational  technology  developed  here. 
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1.  INTRODUCTION 


1.1  Problem  Description 

A  dump  combustor  Is  a  sudden  expansion  which  makes  use  of  the 
sudden  increase  in  area  to  produce  a  flame  holding  effect.  The  flow 
in  dump  combustors  has  been  solved  by  a  number  of  investigators  over 
the  past  few  years.  Examples  of  these  solutions  can  be  found  in 
Gosraan  et  al.  (1),  Lilley  and  Rhode  (2),  Novick  et  al.  (3),  and  Syed 
and  Sturgess  (4).  A  common  characteristic  of  these  solutions  is  the 
mathematical  modeling  of  the  flow  using  the  Navier-Stokes  equations. 

The  flow  in  a  dump  combustor  may  be  viewed  as  having  two  distinct 
regions  with  very  different  characteristics.  The  first  region,  which 
starts  at  the  sudden  expansion  and  extends  to  the  point  of 
reattachraent ,  is  the  recirculation  region.  Because  of  the  presence  of 
the  separation  region  the  pressure  field  exerts  a  strong  upstream 
influence  and  therefore  the  flow  is  said  to  be  elliptic. 

The  second  region,  after  reattachment,  i3  the  pipe  flow  region. 
In  this  region  the  flow  is  "one-way"  (5)  meaning: 

a)  there  is  no  backflow  in  the  main  direction  of  flow, 

b)  the  streamwise  diffusion  of  energy,  momentum,  and  mass  is 
negligible,  and 

c)  the  upstream  flow  conditions  are  unaffected  by  a  downstream 
pressure  field. 

In  thi3  region  the  elliptic  effect  of  the  Navier-Stokes  equations 
disappears  and  the  flow  becomes  parabolic  in  nature.  The  purpose  of 
this  thesis  is  to  describe  a  new  computational  method  for  dump 
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combustor  calculations  which  takes  advantage  of  this  characteristic. 
In  particular  it  is  proposed  to  replace  the  solution  of  the  NavJer- 
Stokes  equations  with  the  solution  of  the  parabolized  Navier-Stokes 
equations  in  this  region.  It  will  be  shown  chat  this  results  in  a 
considerable  savings  of  computational  time  and  computer  storage. 
References  involving  solutions  of  the  parabolized  Navier-Stokes 
equations  are  numerous  and  can  be  found,  for  example,  in  Patankar  and 
Spalding  (5,6)  and  Spalding  (7). 

The  solution  technique  for  the  elliptic  and  parabolized  Navier- 
Stokes  equations  are  very  different.  The  elliptic  method  is  implicit 
so  each  dependent  variable  calculated  in  the  flow  field  will  have  a 
value  stored  at  each  computational  node.  A  line  relaxation  technique 
is  used  to  calculate  the  flow  field  iteratively  until  the  flow  field 
calculated  is  within  a  desired  tolerance.  A  disadvantage  of  using  the 
Navier-Stokes  equations  is  the  large  computational  time  required  to 
carry  out  these  iterations.  Often  a  grid  that  expands  in  the  axial 
directon  is  used  to  decrease  not  only  the  computational  time  but  the 
computer  storage.  Sometimes  the  consequence  of  using  this  coarser  x- 
grid  mesh  is  the  presence  of  a  serious  decrease  in  accuracy. 

The  parabolized  Navier-Stokes  equations  are  obtained  from  the 
Navier-Stokes  equations  by  neglecting  the  streamwise  diffusion 
term3.  This  allows  the  equations  to  be  solved  by  a  streamwise 
marching  integration  procedure.  Clearly  only  a  small  amount  of 
computer  storage  is  needed.  Since  Che  flow  field  is  not  obtained  by 
iteration  the  computation  time  is  3hort. 


The  combined  elliptic  and  parabolic  procedure  just  described  I 
call  the  hybrid  method.  The  purpose  of  this  thesis  will  be  to  discuss 
how  this  hybrid  method  was  developed,  to  discuss  the  interfacing  which 
must  take  place  between  the  two  methods  of  solution,  and  to  compare 
the  results  of  the  hybrid  solution  with  the  elliptic  calculations  and 


experimental  data 


2.  CODE  DEVELOPMENT 


An  elliptic  code  called  STARPIC,  written  by  Ulley  and  Rhode  (2), 
is  used  to  calculate  the  flow  field  in  the  entrance  of  the  dump 
combustor*  This  elliptic  code  solves  the  Navier-Stokes  equations  on  a 
staggered  grid  control  volume  shown  in  Figure  1*  The  TEACH  code  and 
SIMPLE  algorithm  described  in  Patankar  (.8)  are  used*  This  program  was 
made  available  to  me  by  the  Ramjet  Technology  Branch,  Aero  Propulsion 
Laboratory,  at  the  Wright-Patterson  Air  Force  Base*  A  parabolic  code, 
however,  was  not  available,  therefore  a  considerable  amount  of  time 
was  spent  developing  such  a  code  patterned  after  STARPIC*  This 
section  will  be  mostly  concerned  with  the  development  of  this  code* 
The  calculations  are  restricted  to  incompressible,  axisymmetric, 
steady  flow. 

2.1  Differential  Equations 

The  set  of  governing  equations  for  the  parabolic  method  are 
conservation  of  mass,  the  axial  Navier-Stokes  equation  (axial  momentum 
equation),  the  transport  equation  for  turbulent  kinetic  energy,  and 
Che  transport  equation  for  turbulent  dissipation  rate*  These 
equations  are  presented  below  in  cylindrical  coordinates  (5)* 
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Fig.  1.  Interior  Elliptic  Control  Volume 


Continuity: 


i  f—  (rpu)  _3_(rpv)  »  _ 
r  '■3x  3r  > 


(1) 


Axial  Momentum: 


1  rj  (rpuu)  3  (rgvu) 

r  L  3x  3r 


3  f  .u  3u\ 

n  (rr  n> 


(2) 


Turbulent  Kinetic  Energy,  k: 


A  r3  (rpuk)  3  (rpvk) 

r  L3x  3r 


3  ,  k  3k-, 

17  Ct  -)■ 


3r  - 


*slr 

-  r 


(3) 


Turbulent  Dissipation  Energy  Rate,  e: 


A  r3  (rpue)  _3_(rpve)  _3_  r  e  3e> 

r  l3x  3r  ”  3r  Lri  3rJ 


S£ 


(4) 


The  parabolized  computational  procedure  means  that  the  diffusion  terms 


1  3  f  -U  3u  >  1  3  r  nlc  3k  >  1  3  r  _€  36^  /e\ 

717  (cr  *)'  717  (rr  17)'  “d717(rr  17>  (5> 


normally  present  if  the  elliptic  calculations  had  been  carried  out 
have  been  eliminated  from  the  above  equations*  The  rational  for 
eliminating  these  terms  i3  that  they  present  a  mechanism  by  which 
upstream  diffusion  may  take  place  which  is  contrary  to  the  definition 
of  a  parabolic  flow*  The  u,  k,  and  e  superscripts  refer  to  specific 
terms  unique  to  the  axial  momentum,  turbulence  kinetic  energy,  and 
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turbulence  dissipation  rate  equations  respective!/*  The 

ra,  r*.  and  rS  represent  the  diffusion  coefficients  of  axial 
momentum,  turbulence  kinetic  energy,  and  turbulence  dissipation  rate 
respectively*  The  diffusion  coefficients  and  the  linearized  source 
terms  appearing  in  these  equations  will  be  defined  once  the  finite 
difference  form  of  these  equations  has  been  developed* 


2*2  Finite  Difference  Equations 

The  finite  difference  equations  for  the  parabolic  method  use 
control  volumes  based  on  a  semi -staggered  grid  system  illustrated  in 
Figure  2*  The  radial  velocity  (v)  is  calculated  midway  between  nodes 
S  and  P  but  stored  on  node  P*  The  axial  velocity  (u),  like  all  other 
dependent  variables,  is  stored  on  node  P*  The  semi -staggered  grid 
improves  and  simplifies  Che  convection  term  calculation  by  calculating 
Che  velocities  at  Che  center  of  each  control  face  where'  needed*  The 
governing  equations  of  the  parabolic  method  (1-4)  may  be  written  in  a 
form  such  that  $  represents  any  property  which  can  be  convected  or 
diffused  (5),  giving: 

i.  rj  .(purfr)  i_(pvr$)  _  _3_  ,  *  24)1  _  c*  l'6'i 

r  l3x  3r  3r  *«ri  3r J  J  K  1 


From  this  equation  $  may  take  on  the  values  1,  u,  k,  and  e  thereby 
recovering  the  conservation  of  mass,  axial  momentum,  turbulence 
kinetic  energy,  and  turbulence  dissipation  rate  equations 
respectively*  A  general  differential  equation  in  $  is  convenient  for 
presenting  Che  finite  difference  equations* 

The  development  of  the  finite  difference  equations  make 
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assumptions  on  the  change  of  $  between  grid  points  (5). 

a)  In  the  axial  direction  $  is  a  stepwise  property.  $  *  ^  at 
the  downstream  face  of  a  control  volume  up  to,  but  not 
including,  the  upstream  face  of  the  same  control  volume. 

b)  The  upstream  and  downstream  face  of  each  control  volume 
falls  on  the  r9  plane.  The  $  associated  with  an  r9  face  of 
a  control  volume  is  stepwise.  This  means  ^  remains  fixed 
on  the  face  of  the  control  volume  but  suddenly  changes 
to  or  when  exiting  radially  to  an  adjacent  control 
volume. 

c)  The  radial  top  and  bottom  face  of  each  control  volume  is  an 
x9  surface.  The  numerical  value  of  $  convected  on  a  radial 
control  volume  face  will  be  the  arithmetic  mean  of  the 
adjacent  radial  $  values. 

d)  In  the  diffusion  of  $  across  the  same  x9  surface  of  a 
control  volume  a  linear  profile  in  <p  between  adjacent 
radial  $  values  will  be  assumed. 

The  finite  difference  equations  based  on  the  above  assumptions  are 
developed  by  integrating  Equation  7  over  the  control  volume,  giving: 


[l] 


(r  Ar  p  u)p^D  <j>p^D  -  (r  Ar  p  u)p^  $Pj(J 
+  (ra  A*  »  v>n,0  -  Cr3  Ax  p  v)gjD  (^La2_J»i2.) 


(r  Ax  r*)  „  ( ^,t|D-rY— 1-0 )  -  (r  Ax  r*)  n  (■ 

a  a,U  Ai  '  3  s,U 


D  D ' 


+  S^r  Ar  Ax 


See  Table  I  for  one  definition  of  Che  diffusion  coefficients,  and 

a  _ 

source  term  S  .  The  coefficient  (r  Ar  p  u)p  Q  multiplied  by  Q  in 
the  axial  convection  term  must  be  rewritten  to  eliminate  the  unknown 
downstream  influence  on  axial  velocity  (5),  To  eliminate  this  problem 
the  continuity  equation  is  written  in  finite  difference  form  and 
integrated  over  the  control  volume.  The  continuity  equation  is  solved 
for  the  troublesome  coefficient  just  discussed. 


(r  Ar  p  u)_  -  (r  Ar  p  u)  n  -  (r  Ax  p  v)  +  (r  Ax  p  v)  (8) 

r»u  r,u  n  n,u  s  s,u 


Making  this  substitution  and  rearranging  the  terms  in  the 
general  $  finite  difference  equation  produces  the  following  general 
discretization  equations. 


*?,*  -  \u  *  +  As,u  *s  +  30 


TABLE  I 


12 


4,u  "  an,u  +  As,u  ■  r  Ar  ^  s? 


+  (r  Ar  p  a)p>u  -  (rn  Ax  p  v)n>u  +  (rg  Ax  p  v)g>0 


(r  Ax  T?)  n  (r  Ax  p  v)  .. 
4  *  n  n » u  n  a ,  u 

Ti.U  AY  "  2 

a 


(r3  ^  F  VtJ  (C3  **  P  V)3,0 


-  (r  Ar  p  u)p>[]  $p>0  +  “  Ar  Ax  S* 

The  symbol  0  used  as  che  second  subscript  on  the  coefficients  of  $  in 
the  above  equations  satisfies  Che  fact  that  coefficients  are  to  be 
evaluated  using  upstream  (U)  conditions  which  of  course  is 
necessitated  by  the  parabolic  nature  of  the  calculations.  This 
general  equation  calces  on  the  following  specific  form  for  ■)>  ■  1,  u,  k 
and  e  where  S^f  and  are  defined  in  Table  I. 

Let  $  *  l,  then 


(r  Ar  p  u)p  Q  - 

(r  Ar  p  u)p>0  +  Crn 

Ax  p 

v)  .  + 
n,U 

Let  $  *  u, 

then 

£  £ 

^n.u  As,u 

“P 

UM  +■  —  U<,  + 

a“  n  a«  s 

*?,0  T»,U 

•S,u 

Lac  f  m  ic,  Chen 


k 


P 


k 
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Let  $  •  e,  then 


*P 


(13) 


(14) 


The  pressure  term  in  the  axial  equations  will  be  discussed  in  the 
boundary  condition  section  (2*3)*  The  $  coefficients  and  the  B^  terms 
of  the  discretization  equations  are  based  only  on  upstream  influence 
as  indicated  by  the  U  subscript*  The  unknown  <j>'s  in  the 

discretization  equations  are  to  be  solved  at  a  downstream  station  by 
the  tri-diagonal  matrix  (TDM)  algorithm* 


2*3  Boundary  Conditions 

In  the  previous  section  the  finite  difference  equations  were 
developed  for  the  interior  of  the  flow*  This  section  describes 
special  forms  of  these  equations  used  on  the  boundary*  Since  the  flow 
is  axisymmecrlc,  only  the  top  half  of  Che  flow  field  is  considered 
with  the  boundaries  being  the  top  wall  and  the  centerline*  The 
special  forms  for  the  axial  velocity,  turbulence  kinetic  energy,  and 
turbulence  dissipation  rate  will  be  presented  collectively  at  the 
centerline  and  individually  at  the  wall*  The  calculation  of  radial 
velocity  will  be  discussed  in  the  next  section. 
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At  the  centerline  since  the  flow  i3  axisymmetric  the  variables  u, 


k,  and  s  have 

a 

zero  radial 

gradient.  The 

coefficients 

,u  »’* 

As,u’  'S.ir  and 

Ae 

S,U 

belong!. -g  to 

the  centerline 

discretization 

equations  are  set  equal  to  zero  to  satisfy  thi3  boundary  condition. 


In  preparation  for  considering  the  conditions  to  be  applied  along 

the  combustor  wall,  it  is  important  to  keep  in  mind  that  the  axial 

velocity,  turbulence  kinetic  energy,  and  turbulence  dissipation  rate 

equations  apply  only  for  high  Reynolds  numbers.  Near  the  wall, 

however,  the  flow  is  dominated  by  the  laminar  viscosity  (10).  To 

provide  appropriate  modeling  in  the  near  wall  region,  u  and  k 

equations  will  be  modified  by  a  wall  function  and  the  e  equation  will 

force  e  to  be  a  fixed  value.  These  near  wall  equations  will  use  the 

control  volume  shown  In  Figure  3  and  the  wall  function  defined  in 

U  1c  £ 

Appendix  A.  The  coefficients  A^  D,  Ajj  and  A^  Q  associated  with 
the  u,  k,  and  e  discretization  equations  will  be  set  equal  to  zero  so 
that  the  influence  of  the  wall  is  handled  entirely  by  the  source  terms 
S  ,  S  ,  and  S  .  The  source  terms  for  each  near  wall  equation 
supplement  the  source  terms  associated  with  these  equations  for  an 
interior  control  volume. 

The  source  term  needed  to  modify  the  axial  velocity 

discretization  equation  at  the  near  wall  region  follows  (9). 

SU  -  /  t  dV  (15) 

V  w 


where 


;  Y  <  11.63 


.  1/4.  1/2  ,  . 

P  K  Cu  k  <up  '  UN) 

loge  (SY+) 


;  Y  >  11.63 


/T 


1  „ 

U _ 


The  source  term  necessary  to  modify  the  near  wall  turbulent 
kinetic  energy  equation  follows  (9) 


S  -  /  (G  -  Cp  p  e)  dV 
V  r 


The  calculations  of  this  expression  may  be  broken  into 


,  Tw  (uP  “  V 

f  G  dV  -  — - - - —  dV 

V  yp 


/v  Cp  p  e  dV  - 


-  S  Wil>* « 


t  is  calculated  from  Equation  17  and  0+  is  calculated  from 


for  Y  <  11.63 


-  log  (5Y+)  for  Y+  >  11.63 


The  value  of  s,  unlike  value  of  k,  reaches  a  much  higher  value  near 


the  wall  than  anywhere  else  (9).  The  control  volume  for  e  that 
extends  up  to  the  wall  is  very  difficult  to  model  due  to  "shear 
ignorance"  (9)  of  e's  near  wall  profile.  To  eliminate  this 
problem  e  is  given  a  fixed  value  which  is  not  a  function  of  Y+  (9). 
The  initial  sublayer  allows  e  to  be  expressed  as 


c3/V/2 

<  Yo 


(22) 


2.4  Solution  Technique 

The  axial  velocity  can  be  obtained  from  the  axial  momentum 
equation  if  the  axial  pressure  is  ’'nown  from  Equation  13.  The  guessed 
pressure  gradient  used  in  thi3  equation  is  obtained  from  the  upstream 
station.  In  general,  the  estimated  flow  rate  based  on  the  guessed 
pressure  will  not  satisfy  conservation  of  mass.  The  following 
equations  developed  in  Appendix  B  will  be  used  to  correct  pressure 


gradient  and  the  axial  velocity 


is 


j 


j 


z  pi  ui  ri  ^i  “  z  pi  ui  rl 

.  [  >1 1  J— ] 

4  z,  pi  DS,  ri  Ari 

i-1  i 


(25) 


The  *  in  the  above  equations  represents  axial  velocity  based  on  the 

previous  pressure  gradient#  A  numerical  procedure  for  stabilizing  the 
axial  velocity  calculation  i3  presented  in  Appendix  Dl# 

The  discussion  of  radial  velocity  delayed  earlier  will  now  be 

presented#  Radial  velocity  is  calculated  from  Equation  12#  The 

boundary  conditions  for  v  are  zero  at  the  wall  and  centerline#  Since 
the  radial  grid  is  highly  nonuniform,  the  control  volumes  vary  in 

size.  In  order  to  reduce  the  numerical  error  of  the  calculations  the 
radial  velocity  is  calculated  inward  from  the  wall  and  outward  from 
the  centerline  to  the  largest  control  volume#  The  radial  velocity 
calculated  by  the  conservation  of  mass  is  stabilized  by  a  numerical 
procedure  presented  in  Appendix  D2# 

The  turbulence  kinetic  energy  is  calculated  from  Equation  14  and 
then  the  turbulence  dissipation  rate  is  calculated  from  Equation  15# 
Prom  the  value  of  k  and  e  the  turbulent  viscosity  is  calculated  from 
Equation  5  in  preparation  for  taking  the  next  downstream  step# 
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3.  PIPE  FLOW  TEST  CASE 


In  order  to  check  out  the  parabolic  code,  a  simple  test  case  was 
run*  The  test  case  consists  of  laminar  flow  of  air  In  a  circular 
pipe*  This  Is  a  well  known  problem  called  Poiseuille  flow  (11)  for 
which  there  is  an  exact  solution* 

3*1  Problem  Description 

The  data  used  In  this  test  case  was: 

Variable 

Dynamic  Viscosity 
Density 
Diameter 
Temperature 

umax 

Thi3  results  in  a  Reynolds  number  of  2000  which,  since  this  is  a 
laminar  calculation.  Is  necessarily  less  than  the  commonly  accepted 
critical  Reynolds  number  for  pipe  flow  of  2300*  The  results  of  the 
analysis  are  represented  by  Che  following  equations: 


Value 

.1711  X  10  N  s/m2 
1.19  kg/m3 
0.0964  m 
293.0  X 
0.30  m/s 


Ac  the  Inlet  to  the  pipe  the  axial  velocity  profile  given  by 
Equation  27  was  input  along  with  the  value  of  the  pressure  gradient 
given  by  Equation  29.  The  value  of  the  radial  velocity  at  the  inlet 
was  set  equal  to  zero*  From  these  calculations  the  axial  velocity, 
the  radial  velocity,  and  the  axial  pressure  gradient  were 
calculated*  The  axial  mesh  spacing  from  the  calculations  was  0.005 
m*  The  pipe  was  0*6  m  long* 

3*2  Results 

The  initial  velocity  profile  along  with  the  velocity  profiles  of 
every  tenth  station  are  plotted  in  Figure  4.  The  ordinate  represents 
the  radial  grid  spacing  from  the  centerline  to  the  wall.  The  abscissa 
indicates  the  magnitude  of  the  axial  velocity  which  increases  from  0.0 
m/s  at  the  wall  to  u___  *  0*30  m/s  at  Che  centerline.  There  is  no 

lUdX 

change  in  the  axial  velocity  profile  from  the  initial  station  to  the 
final  station  just  as  would  be  expected  from  the  exact  solution. 
Therefore  in  Figure  4  only  a  3ingle  curve  is  plotted  which  in  fact 
represents  Che  initial  velocity  as  well  as  the  velocity  profile  at 
each  of  the  120  axial  stations  at  which  calculations  were  made. 

Although  not  specifically  shown  here,  the  radial  velocities 
remain  effectively  zero.  The  radial  velocity  is  between  4  and  5 
orders  of  magnitude  less  chan  u  with  the  sign  being  alternately 
positive  and  negative  with  each  successive  step. 
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Axial  Velocity  Profile  Calculations  for  Fully  Developed 
Pipe  Flow 


Of  particular  Interest  is  the  pressure  gradient  presented  in 
Table  II.  The  pressure  gradient  calculations  are  presented  at  the 
initial  x  station  and  every  following  tenth  axial  station.  The  value 
of  each  axial  pressure  gradient  is  virtually  the  same  as  the  exact 
solution  (-0.008825). 


TABLE  II 


Axial  Pressure  Gradient  Calculations 
at  Every  Tenth  Axial  Step 
(Step  Size  of  0.005  m) 


ation. 

x(a) 

1 

.00 

10 

.05 

20 

.10 

30 

.15 

40 

.20 

50 

.25 

60 

.30 

70 

.35 

80 

.40 

Axial  Pressure  Gradient 

-.008825 

-.008866 

-.008854 

-.008849 

-.008848 

-.008845 

-.008842 

-.008843 

-.008843 


90 


45 


008841 


4.  INTERFACING 


The  ultimate  utilization  of  the  parabolic  code  is  in  conjunction 
with  an  elliptic  code.  This  necessitates  a  procedure  for 
transitioning  from  the  elliptic  calculations  to  the  parabolic 
calculations.  The  radial  velocity,  turbulence  kinetic  energy, 
turbulence  dissipation  rate,  and  the  dynamic  viscosity  are  calculated 
and  stored  in  exactly  the  same  location  in  the  elliptic  and  parabolic 
codes.  Therefore,  no  special  procedure  will  be  needed  to  transfer 
these  dependent  variables  from  the  elliptic  code  as  initial  conditions 
for  the  parabolic  code.  A  special  procedure,  however,  will  be  needed 
in  conjunction  with  the  axial  velocity.  The  axial  velocity  calculated 
upstream  and  downstream  of  the  parabolic  starting  station  from  the 
elliptic  code  is  averaged  to  obtain  a  value  of  axial  velocity  at  the 
station  where  the  parabolic  calculations  begin: 


A  special  procedure  is  also  necessary  to  obtain  the  initial  pressure 
and  axial  pressure  gradient  for  the  parabolic  code.  Since  the 
elliptic  radial  pressure  is  approximately  the  same  at  each  radial 
location,  the  pressure  at  node  12  (a  point  approximately  halfway 
between  the  wall  and  the  centerline)  is  taken  as  the  starting  pressure 
for  the  parabolic  code.  The  axial  pressure  gradient  is  also  obtained 
at  node  12  by  subtracting  the  pressure  of  the  elliptic  station 
immediately  following  the  elliptic/parabolic  interface  from  the 
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pressure  at  the  interface  station  and  dividing  this  difference  divided 
by  the  axial  distance  between  these  two  stations. 

In  order  to  take  full  advantage  of  the  hybrid  method,  the 
parabolic  code  should  start  as  soon  after  the  reattachment  point  as 
possible,  thereby  saving  the  most  computational  time  and  machine 
storage.  Although  such  a  procedure  is  desirable  from  the  standpoint 
of  saving  computational  time  and  reducing  storage,  the  parabolic 
calculations  starting  too  close  to  the  recirculating  region  run  the 
danger  of  decreasing  accuracy  since  the  parabolizing  assumptions 
neglect  the  streamwise  diffusion  terms.  Consequently  an  Interface 
location  mist  be  sought  which  represents  a  desirable  tradeoff  between 
accuracy,  computational  time,  and  machine  storage. 

One  of  the  original  purposes  for  carrying  out  this  work  was  to 
perform  numerical  experimentation  to  provide  some  guidelines  for 
selecting  an  appropriate  starting  location  for  the  parabolic 
calculations.  However,  the  expanding  grid  used  by  the  elliptic  code 
results  in  such  a  coarse  mesh  spacing  downsteam  of  the  reattachment 
point  that  the  elliptic  calculations  of  the  flow  field  in  this  region 
were  poor.  These  results  made  it  impossible  to  carry  out  a  systematic 
Investigation  of  the  proper  location  for  the  interface  of  the  two 
codes.  More  will  be  said  on  this  in  the  next  section. 


5.  RAMJET  COMBUSTOR  CALCULATIONS 


5#1  Geometry  and  Test  Conditions 

In  order  to  test  the  accuracy  of  the  hybrid  calculations  the  flow 
in  the  dump  combustor  was  calculated  and  compared  with  the  flow 
resulting  from  an  elliptic  calculation  and  with  experimental  data. 
The  dump  combustor  is  shown  in  Figure  5.  The  value  of  R^  was  chosen 
to  be  0.0482  m  and  the  value  of  R2  was  chosen  to  be  0.0254  m.  The 
inlet  velocity  (U  )  prior  to  opposition  was  taken  to  be  25.1  m/s.  Air 
is  the  working  fluid.  The  diameter,  density,  dynamic  viscosity,  and 
temperature  for  the  calculations  are  the  same  as  that  used  in  the  pipe 
flow  test  case  (Section  3). 

5.2  Calculations 

In  a  production-type  application  of  the  hybrid  method,  the 
parabolic  calculations  would  be  initiated  by  a  partial  solution  of  the 
ramjet  combustor  with  the  elliptic  code.  For  the  calculations 
described  here,  however,  the  entire  dump  combustor  geometry  was  solved 
with  STARPIC  to  provide  the  initial  conditions  for  the  parabolic 
calculations.  The  elliptic  code  predicts  reattachment  just  after 
station  35  (x  »  0.18  m)  on  an  11%  expanding  grid  which  extends  from  x 
-  0.0  m  (station  1)  at  the  entrance  to  x  *  0.67  m  (station  48)  at  the 
exit.  The  farthest  upstream  station  at  which  the  parabolic 
calculations  could  begin  was  at  station  35.  This  section  will 
describe  these  calculations  and  compare  the  results  of  the  parabolic 
calculations  carried  out  to  station  40  and  43  with  the  elliptic 


calculations  and  experimental  data.  The  parabolic  calculations  use  a 
constant  x  grid  mesh  spacing  of  0.001  m. 

5.3  Results 

The  solution  of  the  equations  presented  in  Section  4  allows  the 
calculation  of  the  axial  and  radial  velocities,  the  local  static 
pressure,  the  turbulence  kinetic  energy,  and  turbulence  dissipation 
rate  at  each  poinc  in  the  flow  field.  In  the  comparison  presented 
here,  however,  only  the  axial  velocity,  and  turbulent  intensity  are 
shown  since  these  are  the  only  quantities  for  which  experimental  data 
was  available.  The  experimental  data  was  obtained  from  Craig  et  al. 
(12)  and  Stevenson  et  al.  (13).  The  axial  locations  at  which  these 
data  were  taken  correspond  to  station  40  and  43  in  the  present 
calculations  (see  Figure  5).  The  calculations  at  station  40  are  shown 
in  Figures  6  and  7,  The  axial  velocity  and  turbulent  intensity  are 
nondimensionalized  by  the  velocity  just  before  expansion  (UQ)  and 

plotted  as  the  abscissa  in  Figures  6  and  7  respectively.  The  ordinate 
of  both  figures  represents  the  radius  nondimensionalized  by  the  radius 
of  the  pipe  before  expansion  (I^)* 

The  results  in  Figure  6  indicate  as  expected  an  increase  in 
velocity  as  the  centerline  is  approached.  Clearly  the  presence  of  the 
recirculating  zone  is  strongly  felt  at  this  station  and  as  expected  is 
far  from  a  fully  developed  turbulent  profile.  It  is  interesting  that 
Che  parabolic  calculations  are  in  much  better  agreement  with  the 

experimental  data  than  the  elliptic  calculations.  The  result  is  quite 
encouraging.  There  13 ,  however,  a  disappointing  nature  to  this 
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discovery.  It  was  initially  planned  to  use  station  40  to  begin  a 

second  set  of  parabolized  calculations.  The  poor  agreement  between 

the  elliptic  results  and  the  experimental  results  at  station  40 
clearly  point  out  the  futility  of  such  calculations.  Therefore 

parabolic  calculations  beginning  at  station  35  will  be  the  only 

calculations  shown  in  this  thesis.  The  excellent  performance  of  the 
parabolic  calculations  at  station  40  compared  with  the  elliptic 

calculations  requires  an  explanation.  Since  the  elliptic  code  used 
the  full  Navier-Stokes  equations  and  the  parabolic  code  used  an 

abridged  form  of  the  Navier-Stokes  equations,  how  is  it  that  the 
parabolic  code  out  performs  the  elliptic  code?  This  poor  prediction 
of  the  flow  field  from  the  elliptic  code  is  blamed  on  the  highly 
expanded  elliptic  grid.  In  Figure  5A  72.9Z  of  the  nodes  occur  in  the 
first  28.  3Z  of  the  combustor.  This  high  resolution  allows  an  accurate 
description  of  the  flow  field  to  be  obtained  at  station  35.  Between 

station  35  and  48,  however,  the  remaining  71.7Z  of  the  combustor  has 
only  27. 1Z  of  the  nodes.  This  highly  expanding  x  grid  is  blamed  for 
the  severe  loss  of  accuracy  by  the  elliptic  code  in  predicting  the 

flow  field  around  station  40.  In  contrast,  the  success  of  the 
parabolic  calculations  is  attributed  to  the  very  high  axial  grid 
resolution  of  the  calculations. 

In  Figure  7  the  results  of  the  turbulent  intensity  calculations 
are  shown.  Here  even  though  the  agreement  of  the  parabolic 

calculations  of  experimental  data  is  not  as  good  as  the  case  of  the 
axial  velocity  they  still  agree  much  better  with  the  experiments  than 
do  the  elliptic  calculations.  The  poor  agreement  of  the  elliptic 
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calculations  with  the  experimental  turbulent  intensity  again  points 
out  the  futility  of  choosing  station  40  to  begin  the  parabolic 
calculations. 

Figures  8  and  9  are  similar  to  Figures  6  and  7,  however,  the 
calculations  and  experimental  data  are  compared  at  station  43  rather 
than  at  station  40.  The  axial  velocity  and  turbulent  intensity  are 
shown  on  the  same  axes  and  with  the  same  nondimens ional  parameters  as 
in  Figures  6  and  7.  By  the  time  this  station  is  reached  the 
distortions  of  the  flow  produced  by  the  recirculating  zone  are  largely 
diffused  and  dissipated.  The  resulting  flow  calculations  closely 
resemble  turbulent  pipe  flow.  Although  from  Figure  8  it  is  difficult 
to  judge  whether  the  better  performer  is  the  elliptic  or  the  hybrid 
calculation,  in  Figure  9  it  is  quite  clear  that  the  hybrid  results  are 
superior  in  terras  of  predicting  the  turbulent  intensity.  The  results 
which  we  have  just  seen  are  certainly  gratifying  from  the  standpoint 
of  accuracy. 

Given  the  inability  of  the  elliptic  calculation  at  station  40  and 
43  to  adequately  predict  the  flow  field  one  may  be  justifiably 
suspicious  of  the  accuracy  of  the  elliptic  calculations  at  station  35 
where  the  parabolic  calculations  began.  An  examination  of  Figure  10 
should  alleviate  this  concern.  Figure  10,  unlike  Figure  9,  shows  only 
the  turbulent  intensity  on  the  centerline.  As  previously  indicated  at 
stations  40  and  43  the  agreement  is  poor,  however,  at  station  35  where 
the  parabolic  calculations  begin  the  agreement  of  the  experimental 
data  with  the  elliptic  calculations  is  excellent. 
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Although  no  experimental  data  is  available  it  is  nevertheless 
interesting  to  look  at  the  calculations  of  radial  velocity*  The 
results  are  shown  in  Figure  11  where  the  radial  velocity  profile  is 
plotted  as  a  function  of  radius  at  station  35  and  station  43.  In 
addition,  calculations  at  every  twentieth  station  are  plotted  between 
station  35  and  43*  The  results  are  reasonable*  A3  station  43  is 
approached  the  radial  velocity  continually  decreases*  This  is  to  be 
expected  since  the  region  of  fully-developed  pipe  flow  is  being 
approached*  Another  method  for  calculating  the  radial  velocity  can  be 
found  in  the  Appendix*  It  use3  the  radial  momentum  (Havier-Stokes) 
equation*  This  method  is  discussed  in  Appendix  C  with  some 

stabilizing  numerical  procedure  developed  in  Appendix  D3* 

5.4  Computational  Savings 

The  original  motivation  for  the  hybrid  method  was  to  save 
computational  time  and  machine  storage  but  still  maintain  an 
acceptable  degree  of  accuracy*  The  amount  of  computational  time  saved 
will  be  considered  first* 

Let  us  consider  the  requirements  for  the  hybrid  calculations 
first*  The  elliptic  code  which  must  be  run  to  provide  initial 
conditions  to  the  parabolic  calculations  was  run  on  the  grid  shown  in 
Figure  5A*  This  required  3.88  minutes*  The  parabolic  calculations 
required  0.17  minutes*  The  total  time  for  the  hybrid  method  was 
accordingly  4*05  minutes* 

Let  us  now  turn  to  the  time  required  to  make  a  similiar  elliptic 
calculation.  In  order  to  present  a  fair  comparison  it  is  necessary  to 


Radial  Velocity  Plotted  at  Every  Twentieth  Ax  From  Station  35 
to  Station  43  by  the  Hybrid  Method  with  a  Step  Size  of  0.001  m 
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run  our  elliptic  calculations  on  the  grid  shown  in  Figure  5B.  This  in 
fact  was  never  done*  Therefore  an  estimate  for  the  required 
computational  time  was  made*  The  required  computational  time  will  be 
calculated  by  multiplying  the  time  required  to  solve  the  flow  field 
elliptically  using  the  grid  shown  in  Figure  5A  (3.88  minutes)  by  the 
ratio  of  the  number  of  axial  stations  used  in  the  hybrid  method  to  the 
number  of  axial  stations  used  in  Figure  5A  (523/48)*  The  resulting 
estimate  for  the  elliptic  calculations  on  the  same  grid  as  that  used 
by  the  hybrid  method  is  42*3  minutes*  This  estimate  of  computer  time 
is  admittedly  conservative  because  it  is  a  well  known  characteristic 
of  iterative  schemes  that  the  computational  time  does  not  increas  at  a 
linear  rate  with  mesh  points  but  rather  increases  with  the  number  of 
mesh  points  squared  and  in  some  cases  cubed.  As  can  be  clearly  seen 
by  comparing  the  estimated  42*3  minutes  required  by  the  elliptic 
method  to  the  4.05  minutes  required  by  the  hybrid  method  there  is  at 
least  an  order  of  magnitude  savings  in  computational  time  by  using  the 
hybrid  method. 

Let  us  now  consider  the  machine  storage  saved  by  using  the  hybrid 
method  instead  of  the  elliptic  method*  This  is  accomplished  by  first 
calculating  the  machine  storage  of  the  hybrid  method*  The  required 
storage  is  that  of  the  elliptic  code  (needed  as  initial  conditions  for 
the  parabolic  calculations)  and  the  parabolic  calculations 
themselves*  The  machine  storage  used  by  the  elliptic  code  is  obtained 
from  aultiplying  the  5  dependent  variables  (u,  v,  p,  k,  and  e)  by  the 
radial  locations  (24)  and  the  48  axial  planes.  The  elliptic  code 
therefore  requires  5752  storage  locations*  The  machine  storage  needed 
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for  Che  parabolic  code  Is  determined  by  multiplying  Che  4  dependent 
variables  (u,  v,  k,  and  e)  by  the  radial  locations  (24)  and  the  (only) 
2  axial  stations  (an  upstream  and  a  downstream  station).  This  results 
in  192  storage  locations.  Therefore,  a  total  of  5952  machine  storage 
locations  are  needed  for  the  hybrid  method. 

To  compare  with  the  machine  storage  calculated  for  the  hybrid 
method  the  same  grid  system  will  be  used  to  obtain  the  machine  storage 
needed  to  calculate  the  flow  field  with  the  totally  elliptic  method. 
An  estimate  of  the  number  of  storage  locations  is  obtained  by 
multiplying  the  5  dependent  variables  (u,  v,  p,  k,  and  e)  by  the  24 
radial  locations  and  the  523  axial  stations  used  in  the  hybrid  mesh 
shown  in  Figure  5B.  This  multiplication  shows  a  need  for  62,760 
machine  storage  locations  using  the  totally  elliptic  method.  A 
comparison  of  the  5952  machine  storage  locations  needed  for  the  hybrid 
calculations  shows  over  an  order  of  magnitude  savings  in  the  machine 
storage  of  the  hybrid  method. 


6.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  accuracy  at  which  che  paraboLic  code  predicted  the  flow  field 
is  very  satisfying.  Unfortunately,  a  detailed  study  of  the  optimum 
location  of  the  interfacing  point  could  not  be  carried  out;  however,  a 
satisfactory  calculation  was  made.  Savings  in  computational  time  and 
machine  storage  of  more  than  one  order  of  magnitude  was  achieved  with 
no  preceivable  decrease  in  the  accuracy  of  the  calculations  compared 
with  elliptic  calculations  made  with  the  same  axial  resolution. 

By  way  of  a  recommendation  I  feel  that  the  comparisons  just 
reported  should  be  repeated  with  adjustments  made  to  the  elliptic  code 
to  reduce  the  mesh  spacing  downstream  of  the  reattachment  point.  This 
would  allow  the  optimal  starting  location  of  the  parabolic  code  to  be 
found.  The  parabolic  code  should  al3o  be  extended  to  handle 


compressible  and  reacting  flows 
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APPENDIX  A 

The  Wall  Function  Definition 


In  the  near  wall  region  the  local  Reynolds  number  Y+  increases 
rapidly*  In  order  to  accurately  model  the  wall  shear  stress  in  this 
region  it  is  necessary  therefore  to  define  a  wall  function.  The 
development  of  the  wall  function  is  guided  by  the  one-dimensional 
Couette  flow  analysis  in  the  region  close  to  the  wall  where  the  shear 
stress  can  be  assumed  to  be  a  constant.  The  Navier-Stokes  (axial 
momentum)  equation  can  be  reduced  by  these  assumptions  to  the 
following  very  simple  nondimens ional  form  (9). 


T 

T 

W 


*  f)  f* 

l  dy 


(Al) 


where 


and 


Y 


+ 
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(A3) 


The  general  wall  functions  based  on  Y+  are  broken  into  three  regions: 

+  Ui 

a)  Viscous  sublayer  (0  <  Y  <  5)  where  —  »1. 
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b)  Buffer  zone  (5  <  Y+  <  30)  where  laminar  viscosity  and 

turbulent  viscosity  are  of  the  same  magnitude  such  that 
neither  dominates  the  flow. 

u 

c)  Outer  sublayer  (30  <  y+  <  400)  where  —  «  1. 

The  approach  taken  by  STARPIC  and  adopted  by  the  parabolic  scheme  is 
to  eliminate  the  buffer  zone  by  assuming  the  flow  is  purely  viscous 
from  the  wall  to  Y+  <  *  11.63.  For  Y+  >  11.63,  the  flow  is  assumed  to 
be  fully  turbulent  (2).  Y+  undergoes  a  step  change  between  two  wall 
functions  (9).  Equation  A1  may  be  integrated  in  these  two  regions  to 
obtain  a  calculation  of  U+  as  a  function  of  Y+. 
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Y+  ,  if  Y+  <  11.63 

< 

~  loge  (EY+),  if  Y+  >  11.63 


(A4) 


The  integration  constant,  E  ■  9.793,  depends  on  shear  3tress  in  the 
viscous  sublayer  and  pipe  roughness.  The  von  Karman  constant  <  * 
0.4187. 

In  the  inertial  region  (0  <  Y+  <  400)  mean  velocity  is  considered 
to  be  independent  of  7  so  the  time  averaged  convection  and  diffusion 
terms  in  the  energy  transport  equations  average  to  zero.  The 
turbulence  equations  reduce  to  Production  -  Dissipation  (14).  This 
assumption  in  the  inertial  region  formulates  the  basis  for  the  near 
wall  shear  stress  (t  )  calculations  as  a  function  of  turbulent  kinetic 
energy.  An  isotropic  viscosity  assumption  (2)  further  reduces  the  k 


equation  to 
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US 

Correction  for  Axial  Velocity  and  Prassure 

A  pressure  correction  equation  is  needed  to  correct  the  axial 
velocity  and  update  the  pressure  gradient.  Equation  2,  rewritten 
below  for  convenience,  is  used  to  calculate  the  axial  velocity  with 
the  correct  pressure  gradient. 


*£,0  AS 


However,  the  pressure  gradient  is  unknown  so  a  guessed  pressure 
gradient  will  be  used  to  obtain  an  estimate  to  the  axial  velocity  as 
indicated  by  the  *  in  the  following  equation. 

u  u  u  u  _  * 

*  Si  o  *  AS.U  *  BU  u 


In  general  the  axial  velocity  profile  obtained  from  this  guessed 
pressure  field  will  not  satisfy  the  following  conservation  of  mass 
equation. 


2u  Z  p.  u,  r,  Ar. 
i-1  1  1  L 


The  difference  between  Equation  31  and  Equation  B2  would  represent  the 
error  associated  with  calculating  the  axial  velocity  from  the  guessed 


prassure  as  shown  by 
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*  *  Vu  r  * >  .  Xu  r 

“p  +  C«s  -  v  +  7^“  ius 
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where 


3P 

3x 


(B5) 


For  the  above  equation  to  functioa  as  a  velocity  correction  equation 
it  oust  be  explicit,  therefore,  the  second  and  third  terms  on  the 
right  are  eliminated. 
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The  explicit  axial  velocity  correction  equation  shown  above  is 
substituted  into  Equation  B3  to  calculate  the  following  pressure 
gradient  correction. 
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This  equation  gives  the  required  (— }  necessary  to  update  the 

3P*  * 

guessed  pressure  gradient  (— )  from  Equation  B5  and  to  correct  the 
axial  velocity  from  Equation  B6. 
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APPENDIX  G 


Radial  Momentum  Equation  for  v 


An  alternative  to  the  conservation  of  mass  equation  for  the 
calculation  of  the  radial  veocity  (v)  is  the  following  radial  momentum 
(Navier-Stokes )  equation: 


1  r 3_( r puv )  3  (rpw)  3  r  v  3vn  i 
r  L3x  3r  ”  3r  3v'J 


(Cl) 


The  radial  momentum  equation  may  be  written  in  finite  difference  form 
and  integrated  over  a  control  volume  to  produce  the  following 
discretization  form: 


As  in  the  calculation  of  v  from  the  conservation  of  mass  the  radial 
velocity  calculations  follow  the  axial  velocity  and  pressure  gradient 
calculations*  Since  the  radial  pressure  field  is  unknown  a  guessed 
radial  pressure  field  is  substituted  into  Equation  C2  to  obtain  an 
estimate  for  v  as  shown  below* 


(C3) 


The  *  indicates  the  calculated  radial  velocity  is  based  on  a  guessed 
radial  pressure  field  (?*)•  The  radial  velocity  based  on  the  P*  will 
not  satisfy  in  general  the  conservation  of  mass  equation  by  producing 
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a  net  gain  or  loss  in  mass*  A  procedure  Co  annihilate  this  net  mass 
is  to  create  a  pressure  correction  from 


P 


* 


P  +  P 


(C4) 


I 

where  P  is  the  pressure  correction,  P  is  the  guessed  pressure,  and  P 
is  the  actual  pressure*  The  radial  velocity  correction  equation  will 
be  derived  by  subtracting  the  actual  radial  velocity  equation 
(Equation  G2)  based  on  the  actual  radial  pressure  field  from  the 
radial  velocity  equation  (Equation  C3)  based  on  the  guessed  radial 
pressure  field*  The  following  radial  velocity  correction 
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is  reduced  to 
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(C6) 


by  dropping  the  last  two  terms*  The  last  two  terms  must  be  dropped  in 
order  to  reduce  Equation  C5  from  enplicit  to  explicit  Equation  C6. 
The  above  corrected  radial  velocity  equation  is  substituted  into  the 
finite  differenced  conservation  of  mass  equation  (Equation  12)  to 

I 

develop  the  following  P  correction  equation 
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This  aquation  is  solved  by  a  TDM  procedure  to  calculate  the  P',s« 
With  the  P'.s  known,  the  radial  pressure  may  be  updated  from  Equation 
C4  and  the  radial  velocity  may  be  corrected  from  Equation  C6. 

The  dump  combustor  calculations  described  in  Section  5  were 
recalculated  with  the  radial  velocity  obtained  from  the  radial 
momentum  equation  instead  of  the  conservation  of  mass  equation*  The 
radial  velocity  calculation  were  essentially  the  same  as  that 
calculated  from  the  conservation  of  mass* 
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Some  Stabilizing  Numerical  Procedures 

The  calculations  of  axial  velocity  and  radial  velocity  required 
some  stabilizing  numerical  procedures  suggested  by  Mrs*  Joan  Moore* 


1*  Axial  Momentum 

The  calculation  of  axial  velocity  obtained  from  the  solution  tech¬ 
nique  section  was  found  to  oscillate  about  the  correct  axial  velocity 
from  one  axial  step  to  the  next*  The  axial  velocity  oscillation  was 
eliminated  by  dividing  the  pressure  correction  term  in  Equation  25  by  2 
as  shown  below: 


UP 


2,  Continuity 

1)  When  the  expressions 
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and 
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are  substituted  into  the  following  finite  differenced  conser¬ 
vation  of  mass  equation 
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and  solved  for  vD,  a  numerical  instability  results.  This 
problem  is  eliminated  by  assuming  (pv)^  and  (pv)  are  located 
at  downstream  stations. 

2.  The  radial  velocity  was  obtained  from  a  highly  nonuniform 
radial  grid  by  solving  the  conservation  of  mass  equation  for 
radial  velocity  from  the  wall  to  the  centerline.  At  the 
centerline  the  numerical  round  off  error  was  significant.  The 
solution  was  to  solve  the  conservation  of  mass  equation  for 
radial  velocity  by  calculating  inward  from  the  wall  and  outward 
from  the  centerline  to  the  largest  control  volume.  The  largest 
control  volume  absorbed  the  numerical  error. 

3.  The  radial  velocity  was  found  to  oscillate  about  the  correct 
radial  velocity  profile  form  one  axial  step  to  the  next.  This 
radial  velocity  oscillation  was  eliminated  by: 


(D5) 


3.  Radial  Momentum 

l.  A  pressure  correction  equation  needed  to  correct  the  radial 
velocity  is  obtained  for  each  control  volume.  However,  since  a 
global  pressure  correction  has  already  been  obtained  to  correct 


the  axial  velocity  and  update  the  axial  pressure  a  full  set  of 
pressure  correction  equations  would  be  indeterminate*  This  is 
due  to  the  presence  of  one  too  many  pressure  correction 
equations.  To  eliminate  this  problem  the  pressure  correction 
associated  with  the  largest  control  volume  is  set  equal  to 
zero* 

The  radial  velocity  was  found  to  oscillate  about  the  correct 
radial  velocity  at  each  axial  step*  The  correction  is  the  same 
as  that  developed  earlier  in  Equation  D5* 
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